{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Isotropic linear-elastic media\n",
    "\n",
    "In the last chapter we derived the equations of motion for wave propagation in a general anisotropic linear-elastic medium. Here, we will derive the isotropic linear-elastic case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stress-Strain relation for isotropic linear-elastic media\n",
    "\n",
    "The equations of motion for a general anisotropic linear-elastic medium consist of the conservation of momentum  \n",
    "\n",
    "\\begin{equation}\n",
    "\\rho \\frac{\\partial^2 u_i}{\\partial t^2} = \\sum_{j=1}^{3} \\frac{\\partial}{\\partial x_j} \\sigma_{ij} + f_i\n",
    "\\end{equation}\n",
    "\n",
    "and the stress-strain relationship\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\sum_{k=1}^{3}\\sum_{l=1}^{3} c_{ijkl}\\; \\epsilon_{kl}\n",
    "\\end{equation}\n",
    "\n",
    "In the snapshots below, we can compare the pressure wavefield in a homogeneous, anisotropic medium with an isotropic medium\n",
    "\n",
    "<img src=\"images/ISO_vs_ANISO.png\" style=\"width: 800px;\"/>\n",
    "\n",
    "In an anisotropic medium, the seismic velocities depend on the direction in which the wave is propagating. From the deformation of the seismic wavefront, you can derive that the P-wave velocity in horizontal direction seem to be larger than in vertical direction. \n",
    "\n",
    "For comparison, in the isotropic elastic medium, the wavefront has a rotational symmetry, because the elastic parameters are not directional dependent.\n",
    "\n",
    "To introduce isotropy we only have to modify the elastic tensor $c_{ijkl}$, so eq.(1) is unaffected. In an isotropic body, which means uniformity in all directions, only 2 elastic constants are independent and the stiffness tensor $c_{ijkl}$ takes the form (Hodge 1961, Jeffreys 1969):\n",
    "\n",
    "\\begin{equation}\n",
    "c_{ijkl}=\\delta_{ij} \\delta_{kl} \\lambda + (\\delta_{ik} \\delta_{jl} + \\delta_{il} \\delta_{jk}) \\mu\n",
    "\\end{equation}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda\\; \\text{and}\\; \\mu:\\; &\\text{Lamé parameters} \\nonumber \\\\ \n",
    "\\delta_{ij}=\n",
    "\\begin{cases}\n",
    "1 & \\rm{\\text{if } i = j}\\\\\n",
    "0 & \\rm{\\text{if } i \\ne j}\n",
    "\\end{cases} : \\;&\\text{Kronecker's delta} \\nonumber\n",
    "\\end{align}\n",
    "\n",
    "Mathematically, eq. (3) states that the components of the elastic tensor remain the same under a rotation of the coordinate system.\n",
    "\n",
    "To calculate the stress-strain relation for an isotropic linear-elastic medium, we have to insert eq. (3) into (2). \n",
    "\n",
    "##### Exercise \n",
    "For given values of i and j, you can compute the components of the stress tensor manually, e.g. $\\sigma_{13}$ by setting i=1 and j=3 and evaluate the sums and Kronecker Deltas in eqs. (3) & (2). However, this is a quite time consuming process.\n",
    "\n",
    "Alternatively, we can use the symbolic math library **SymPy**. To install SymPy you have to open a terminal and type\n",
    "\n",
    "``conda install sympy`` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import SymPy\n",
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test symbolic KroneckerDelta function\n",
    "KroneckerDelta(1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimate stress-strain relation for isotropic linear-elastic medium \n",
    "# using given indices i,j\n",
    "def stress_strain_iso(i,j):\n",
    "    \n",
    "    # shift indices by -1\n",
    "    i-=1\n",
    "    j-=1\n",
    "    \n",
    "    # define cijkl for isotropic medium\n",
    "    c_1 = KroneckerDelta(i,j)*KroneckerDelta(k,l)\n",
    "    c_2 = KroneckerDelta(i,k)*KroneckerDelta(j,l)\n",
    "    c_3 = KroneckerDelta(i,l)*KroneckerDelta(j,k)\n",
    "    cijkl = c_1 * lam + (c_2 + c_3) * mu\n",
    "\n",
    "    # sum over l-index\n",
    "    sum_l = Sum(cijkl*e[k,l], (l, 0, 2)).doit()\n",
    "\n",
    "    # sum over k-index to estimate stress\n",
    "    stress = Sum(sum_l, (k, 0, 2)).doit()\n",
    "    \n",
    "    return stress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define indices\n",
    "i, j, k, l = symbols('i j k l', integer=True)\n",
    "\n",
    "# define components of the deformation tensor as symbols\n",
    "e11, e12, e13, e22, e23, e33 = symbols('e11 e12 e13 e22 e23 e33')\n",
    "\n",
    "# define material parameters\n",
    "lam, mu = symbols('lambda mu')\n",
    "\n",
    "# define deformation tensor\n",
    "e = Matrix([[e11,e12,e13], [e12,e22,e23], [e13,e23,e33]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 e_{13} \\mu\n"
     ]
    }
   ],
   "source": [
    "# stress-strain relation\n",
    "stress = stress_strain_iso(1,3)\n",
    "print(latex(stress))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This was a nice introduction to the computing capabilities of **SymPy**. However, instead of rushing to evaluate the stress-strain relation for an isotropic medium by explicitly calculating the sums over the Kronecker Deltas for each stress component, let's take a closer look again, what we have to evaluate:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\sum_{k=1}^{3}\\sum_{l=1}^{3} \\{\\delta_{ij} \\delta_{kl} \\lambda + (\\delta_{ik} \\delta_{jl} + \\delta_{il} \\delta_{jk}) \\mu\\} \\epsilon_{kl} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "First, we rewrite the relation\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\sum_{k=1}^{3}\\biggl\\{\\sum_{l=1}^{3} \\delta_{ij} \\delta_{kl} \\epsilon_{kl} \\lambda + \\mu \\sum_{l=1}^{3} \\delta_{ik} \\delta_{jl} \\epsilon_{kl} + \\mu \\sum_{l=1}^{3} \\delta_{il} \\delta_{jk} \\epsilon_{kl} \\biggr\\}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and evaluate the sums over all l-indices. The term $\\delta_{ij} \\delta_{kl} \\epsilon_{kl} \\lambda$ is only non-zero if $l=k$, so we get:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\sum_{k=1}^{3}\\biggl\\{\\delta_{ij} \\epsilon_{kk} \\lambda + \\mu \\sum_{l=1}^{3} \\delta_{ik} \\delta_{jl} \\epsilon_{kl} + \\mu \\sum_{l=1}^{3} \\delta_{il} \\delta_{jk} \\epsilon_{kl} \\biggr\\}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "The term $\\delta_{ik} \\delta_{jl} \\epsilon_{kl}$ is only non-zero if $l=j$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\sum_{k=1}^{3}\\biggl\\{\\delta_{ij} \\epsilon_{kk} \\lambda + \\mu \\delta_{ik} \\epsilon_{kj} + \\mu \\sum_{l=1}^{3} \\delta_{il} \\delta_{jk} \\epsilon_{kl} \\biggr\\}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "The opposite is true for the term $\\delta_{il} \\delta_{jk} \\epsilon_{kl}$, where only $l=i$ contributes to the summation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\sum_{k=1}^{3}\\biggl\\{\\delta_{ij} \\epsilon_{kk} \\lambda + \\mu \\delta_{ik} \\epsilon_{kj} + \\mu \\delta_{jk} \\epsilon_{ki} \\biggr\\}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Next, we can evaluate the summation over the k-index, where the term $\\mu \\delta_{ik} \\epsilon_{kj}$ is only non-zero for $k=i$ and $\\mu \\delta_{jk} \\epsilon_{ki}$ only non-zero for $k=j$, leading to \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\delta_{ij} \\lambda \\biggl\\{\\sum_{k=1}^{3} \\epsilon_{kk} \\biggr\\} + \\mu  \\epsilon_{ij} + \\mu \\epsilon_{ji} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Due to the symmetry of the strain-tensor $\\epsilon_{ij} = \\epsilon_{ji}$, we get\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij} = \\delta_{ij} \\lambda \\biggl\\{\\sum_{k=1}^{3} \\epsilon_{kk} \\biggr\\} + 2 \\mu \\epsilon_{ij} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "And finally, the general **stress-strain relation for an isotropic linear-elastic medium** as \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{ij}=\\lambda \\theta \\delta_{ij} + 2 \\mu \\epsilon_{ij} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "with the cubic dilatation\n",
    "\n",
    "\\begin{equation}\n",
    "\\theta = tr(\\epsilon_{ij}) =  \\epsilon_{11} + \\epsilon_{22} + \\epsilon_{33} = \\frac{\\partial u_1}{\\partial x_1} + \\frac{\\partial u_2}{\\partial x_2} + \\frac{\\partial u_3}{\\partial x_3} \\nonumber\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D equations of motion for isotropic linear-elastic medium\n",
    "\n",
    "Writing out the momentum conservation and stress-strain relation for the isotropic linear-elastic medium, we get\n",
    "\n",
    "\\begin{align}\n",
    "\\rho\\frac{\\partial^2 u_1}{\\partial t^2} &= \\frac{\\partial \\sigma_{11}}{\\partial x_1} + \\frac{\\partial \\sigma_{12}}{\\partial x_2} + \\frac{\\partial \\sigma_{13}}{\\partial x_3} + f_1 \\nonumber \\\\\n",
    "\\rho\\frac{\\partial^2 u_2}{\\partial t^2} &= \\frac{\\partial \\sigma_{21}}{\\partial x_1} + \\frac{\\partial \\sigma_{22}}{\\partial x_2} + \\frac{\\partial \\sigma_{23}}{\\partial x_3} + f_2 \\nonumber \\\\\n",
    "\\rho\\frac{\\partial^2 u_3}{\\partial t^2} &= \\frac{\\partial \\sigma_{31}}{\\partial x_1} + \\frac{\\partial \\sigma_{32}}{\\partial x_2} + \\frac{\\partial \\sigma_{33}}{\\partial x_3} + f_3 \\nonumber \\\\\n",
    "\\sigma_{11} &= \\lambda(\\epsilon_{11}+\\epsilon_{22}+\\epsilon_{33}) + 2 \\mu \\epsilon_{11}\\nonumber \\\\\n",
    "\\sigma_{22} &= \\lambda(\\epsilon_{11}+\\epsilon_{22}+\\epsilon_{33}) + 2 \\mu \\epsilon_{22}\\nonumber \\\\\n",
    "\\sigma_{33} &= \\lambda(\\epsilon_{11}+\\epsilon_{22}+\\epsilon_{33}) + 2 \\mu \\epsilon_{33}\\nonumber \\\\\n",
    "\\sigma_{12} &= 2 \\mu \\epsilon_{12}\\nonumber \\\\\n",
    "\\sigma_{13} &= 2 \\mu \\epsilon_{13}\\nonumber \\\\\n",
    "\\sigma_{23} &= 2 \\mu \\epsilon_{23}\\nonumber \\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We learned:\n",
    "\n",
    "* 3 different approaches to derive the stress-strain relation for an isotropic linear-elastic medium\n",
    "\n",
    "  1. Manually, by straightforward computation of the stress tensor components\n",
    "  2. Using the symbolic math library **SymPy**\n",
    "  3. General evaluation of the stress-strain relation and neglecting all non-contributing terms\n",
    "* 3D equations of motion for isotropic linear-elastic medium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Hodge, P.G. (1961) \"[On Isotropic Cartesian Tensors](https://www.jstor.org/stable/2311997?seq=1#page_scan_tab_contents)\", The American Mathematical Monthly, 68(8), 793-95.\n",
    "\n",
    "2. Jeffreys, H. (1969) \"[Cartesian Tensors](https://archive.org/details/CartesianTensors)\", Cambridge University Press\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
